# Dickstein, Ho, and Mark (2023)
# This script collects the specification details, estimates demand, and saves it. 

#################
# PRELIMINARIES #
#################
library(tidyr)
library(readr)
library(dplyr)
library(data.table)
library(knitr)
library(ggplot2)
library(readxl)
library(numDeriv)
library(stringr)
library(yaml)
library(gtools)
library(haven)
library(nloptr)

wd <- paste0(specs_location, "/../..")

##################################
#-----**   HOUSEKEEPING  **------#
##################################
load_internal_scripts <- function(wd) {
  source(paste0(wd, "/get_llh.R"))
  source(paste0(wd, "/parametrize_choice_data.R"))
  source(paste0(wd, "/parametrize_cost_data.R"))
}
load_internal_scripts(wd)

# Start log
sink(file = paste0(output_path, "/log.txt"))

# Adjust data for type (This script creates switcher variables if needed and spits out W1_names, W2_names, and Z_names
source(paste0(wd, "/adjustfortype.R"), local = T)

# Reformat data
data <- as.data.frame(data)

# Adjust for high acg
data$sum_concurrent_risk <- ifelse(data$sum_concurrent_risk > sumrisk_cutoff, sumrisk_cutoff, data$sum_concurrent_risk)

##################################
#-----** ESTIMATION CODE **------#
##################################
print(paste0("Censoring cost at ", fixed_cost_censor))

covar_names <- c(W1_names, W2_names, W3_names, Z_names)
covar_lens <- c(length(W1_names), length(W2_names), length(W3_names), length(Z_names))

cost_covar_names <- c(W1_names, W2_names)
cost_covar_lens <- c(length(W1_names), length(W2_names))

# Drop variables that we don't use in estimation
necessary_vars <- c("subscriberid_year", "choice", "p", "cost", "x_av", "insurance")
est_data <- data %>% select(necessary_vars, covar_names)

lb_w1 <- rep(-3, length(W1_names))
lb_w2 <- rep(-20, length(W2_names))
lb_w3 <- rep(-20, length(W3_names))
lb_z  <- rep(-Inf, length(Z_names))
lb <- c(lb_w1, lb_w2, lb_w3, lb_z)

ub_w1 <- rep(5, length(W1_names))
ub_w2 <- rep(5, length(W2_names))
ub_w3 <- rep(0, length(W3_names))
ub_z  <- rep(Inf, length(Z_names))
ub <- c(ub_w1, ub_w2, ub_w3, ub_z)

global_opts <- list(
  "algorithm" = "NLOPT_LD_LBFGS",
  "xtol_rel" = 1.0e-4,
  "check_derivatives" = check_derivatives,
  "check_derivatives_print" = "all",
  "maxeval" = 2000,
  "print_level" = 3)

if(est_type == "full_only" & (2014 %in% yrvec & 2015 %in% yrvec & 2016 %in% yrvec)){
starting <- c(
  -0.243604, -0.427684, -0.292689, -0.653091, 
  -1.325654,
  -2.214164,
  -4.474871, -2.464714, -1.935607, -5.238182, -3.486299, -3.164166, -2.654556, 0.885239)
} else {
  starting <- rep(-.5, length(covar_names))
}

print("Starting values:")
print(starting)

print("Covar names and lengths")
print(covar_names)
print(covar_lens)

print("Upper and lower bounds")
print(ub)
print(lb)

out_0 <- nloptr(x0=starting,
                eval_f=get_combined_llh,
                opts=global_opts,
                lb=lb,
                ub=ub,
                data=est_data,
                covar_names=covar_names,
                covar_lens=covar_lens,
                cost_covar_names=cost_covar_names,
                cost_covar_lens=cost_covar_lens,
                gradient=TRUE,
                scores=FALSE,
                fixed_cost_censor=fixed_cost_censor)
print(out_0)

##################################
#---------** REPORTING **--------#
##################################

est_names = c(
  "out_0"
)
out = list(
  out_0
)
out <- setNames(as.list(out), est_names)

out_solution <- data.frame()
for (name in names(out)) {
  out_solution <- rbind(out_solution, out[[name]]$solution)
}
names(out_solution) <- covar_names
out_solution <- cbind(est_names, out_solution)

save(out, file=paste0(output_path, "/estimation_outputs.rda"))
write.csv(out_solution, paste0(output_path, "/solution.csv"))
write.csv(covar_names, paste0(output_path, "/covar_names.csv"))
write.csv(covar_lens, paste0(output_path, "/covar_lens.csv"))

rm(est_data)
